
#include <iostream>
#include<exception>
#include<eigen3/Eigen/Core>
#include<eigen3/Eigen/Dense>
using namespace std ;

// with schmidt or....

// TODO what about  rectangular matrix? 
void QRDecomposition(Eigen::MatrixXd input,Eigen::MatrixXd& Q,Eigen::MatrixXd& R)
{
    if(input.rows()<input.cols())
        throw runtime_error("please keep rows>= cols");
    
    for (int i = 0 ;i<input.cols() ;++i)
    {
        Q.col(i)=input.col(i); 
        for ( int j =0;j<i;++j)
        {
            double k=input.col(i).dot(Q.col(j))/(Q.col(j).dot(Q.col(j)));
            Q.col(i)-=k*Q.col(j);
        }
        Q.col(i)=Q.col(i)/(sqrt(Q.col(i).dot(Q.col(i))));
    }
    R=Q.inverse() *input;
}


double findMax(Eigen::MatrixXd x)
{
    if (x.rows()==0)
        throw runtime_error("can not input a zero dim vector");
    double max=0;
    for (int i =0 ;i<x.rows();++i)
    {
        if (x(i) >max)
        max=x(i);
    }
    return max;
}


double  powerMethod(Eigen::MatrixXd A,double err)
{
    if(A.rows()!=A.cols())
        throw runtime_error("A is not square matrix !");
    Eigen::VectorXd x(A.rows());
    x.setRandom();
    double max=findMax(x);
    Eigen::VectorXd x1=A*(x/max);
    double max1=findMax(x1);
    while(abs(max1-max)>err)
    {
        x=x1;
        max=max1;
        x1=A*(x/max);
        max1=findMax(x1);
    }
    return max1;
}

int main()
{
    Eigen::MatrixXd A(100,100);
    // A<<4,5,6,7,8,9,13,12,11,12,13,19,21,20,6,7; 
    A.setRandom(); 
    Eigen::MatrixXd Q(100,100),R(100,100);
    QRDecomposition(A,Q,R);
    cout << Q<<endl; 
    cout <<"R"<<endl;
    cout << R;
    return 0;
}